library(raster) # Raster Data handling
library(tidyverse) # Data Manipulation
library(getlandsat) # keyless Landsat data (2013-2017)
library(sf) # Vector data processing
library(mapview) # Rapid Interactive visualization

#Question 1

bb <- read_csv("~/github/geog-176A-labs/data/uscities.csv") %>%
  filter(city == "Palo") %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>% 
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()
plot(bb)

#Question 2

bbwgs = bb %>% st_transform(4326)
bb = st_bbox(bbwgs)

meta = read_csv("~/github/geog-176A-labs/data/palo-flood.csv")

files = lsat_scene_files(meta$download_url) %>% 
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file)) %>% 
  arrange(file) %>% 
  pull(file)

st = sapply(files, lsat_image)
s = stack(st) %>%
  setNames(paste0("band", 1:6))
    plot(s)

s
## class      : RasterStack 
## dimensions : 7811, 7681, 59996291, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 518085, 748515, 4506585, 4740915  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 
## names      : band1, band2, band3, band4, band5, band6 
## min values :     0,     0,     0,     0,     0,     0 
## max values : 65535, 65535, 65535, 65535, 65535, 65535

dimensions : 7811, 7681, 59996291, 6 (nrow, ncol, ncell, nlayers)

resolution : 30, 30 (x, y)

crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

cropper = bbwgs %>%
  st_transform(crs(s))

r = crop(s, cropper)
r = r %>% 
  setNames(c("Coastal Aerosol", "Blue", "Green", "Red", "Near Infrared", "SWIR 1"))

plot(r)

r
## class      : RasterBrick 
## dimensions : 340, 346, 117640, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594075, 604455, 4652535, 4662735  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : Coastal.Aerosol,  Blue, Green,   Red, Near.Infrared, SWIR.1 
## min values :            8566,  7721,  6670,  6057,          5141,   4966 
## max values :           18482, 19079, 19592, 21965,         24772,  23896

dimensions : 340, 346, 117640, 6 (nrow, ncol, ncell, nlayers)

resolution : 30, 30 (x, y)

crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

#Question 3

R-G-B (Natural Color Image)

plotRGB(r, r = 4, g = 3, b = 2)

NIR-R-G(fa) (Traditional Color Infrared (CIR) Image)

plotRGB(r, r = 5, g = 4, b = 3)

NIR-SWIR1-R (False Color good for picking out land from water)

plotRGB(r, r = 5, g = 6, b = 4)

SWIR2-SWIR1-R (False Color useful for visualizing urban environments)

plotRGB(r, r = 7, g = 6, b = 4)

Colored stretch is an effective and simple way to keep a visual track of different products because colored stretch is easier for people to quickly identify and separate different elements, ensuring a specific element isn’t mixed up with others.

RGB image comparison w stretch=“hist”

par(mfrow = c(2,2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "hist")
plotRGB(r, r = 5, g = 4, b = 3, stretch = "hist")
plotRGB(r, r = 5, g = 6, b = 4, stretch = "hist")
plotRGB(r, r = 7, g = 6, b = 4, stretch = "hist")

RGB image comparison w stretch=“lin”

par(mfrow = c(2,2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3, stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4, stretch = "lin")
plotRGB(r, r = 7, g = 6, b = 4, stretch = "lin")

#Question 4

ndvi = (r$Near.Infrared - r$Red) / (r$Near.Infrared + r$Red)
ndwi = (r$Green - r$Near.Infrared) / (r$Green + r$Near.Infrared)
mndwi = (r$Green - r$SWIR.1) / (r$Green + r$SWIR.1)
wri = (r$Green + r$Red) / (r$Near.Infrared + r$SWIR.1)
swi = 1 / (sqrt(r$Blue - r$SWIR.1))
stack = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue","white","red"))
plot(stack, col = palette(256))

All 5 images emphasize the water feature, and all of them are in blue and red. However, the water features are presented in different color. In NDVI and SWI, the threshold of water features are less than 0. In NDWI, MNDWI, and WIR, the threshold of water features are greater than 0. The range of the water threshold is the greatest in MNDWI and the least in SWI

ndvi2 = function(x){ifelse(x <= 0,1, 0)}
ndwi2 = function(x){ifelse(x >= 0,1, 0)}
mndwi2 = function(x){ifelse(x >= 0,1, 0)}
wri2 = function(x){ifelse(x <= 1,1, 0)}
swi2 = function(x){ifelse(x<=5, 1, 0)}

ndvi3 = calc(ndvi, ndvi2)
ndwi3 = calc(ndwi, ndwi2)
mndwi3 = calc(mndwi, mndwi2)
wri3 = calc(wri, wri2)
swi3 = calc(swi, swi2)

threshold = stack(ndvi3, ndwi3, mndwi3, wri3, swi3) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI")) 
palette = colorRampPalette(c("blue","white"))
plot(threshold, col = palette(256))

#Question 5

set.seed(09072020)
values = getValues(r) %>%
  na.omit(values)
dim(r)
## [1] 340 346   6
threshold <- na.omit(threshold)

There are 6 layers with 59996291 values in each layers.

k12 = kmeans(values, 12, iter.max = 100) kmeans_raster = threshold\(ndvi values(kmeans_raster) = k12\)cluster plot(kmeans_raster)

kmeans_vals = getValues(kmeans_raster) binary_vals = getValues(flood_ndvi) table <- table(binary_vals,kmeans_vals)

which.max(table)

t_kmeans = function(x){ifelse(x == 3, 1, 0)} f_kmeans = calc(kmeans_raster, t_kmeans)

threshold = addLayer(threshold, f_kmeans) %>% setNames(c(“ndvi”, “ndwi”, “mndwi”, “wri”, “swi”, “kmeans”)) palette = colorRampPalette(c(“blue”,“white”)) plot(threshold, col = palette(256))

#Question 6

kabletable = cellStats(threshold, sum)
knitr::kable(kabletable, caption = "Number of Flooded Cells", col.names = c("Number"))
Number of Flooded Cells
Number
NDVI 6666
NDWI 7212
MNDWI 11939
WRI 109172
SWI 15201 
areakable = kabletable * 900
knitr::kable(areakable, caption = "Area of Flooded Cells (m^2)", col.names = c("Area"))
Area of Flooded Cells (m^2)
Area
NDVI 5999400
NDWI 6490800
MNDWI 10745100
WRI 98254800
SWI 13680900

Extra Creedit

point = st_point(c(-91.78948,42.06306)) %>% 
  st_sfc(crs = 4326) %>% 
  st_transform(crs(threshold)) %>%
  as_Spatial()
print("-91.78948,42.06306")
## [1] "-91.78948,42.06306"
raster::extract(threshold, point)
##      NDVI NDWI MNDWI WRI SWI
## [1,]    1    1     1   0   1